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Abstract The complexity of the interactions between the constituent gran¬ 
ular and liquid phases of a suspension requires an adequate treatment of the 
constituents themselves. A promising way for numerical simulations of such 
systems is given by hybrid computational frameworks. This is naturally done, 
when the Lagrangian description of particle dynamics of the granular phase 
finds a correspondence in the fluid description. In this work we employ ex¬ 
tensions of the Lattice-Boltzmann Method for non-Newtonian rheology, free 
surfaces, and moving boundaries. The models allows for a full coupling of the 
phases, but in a simplified way. An experimental validation is given by an 
example of gravity driven flow of a particle suspension. 

Keywords Suspensions • Lattice-Boltzmann Method • Discrete Element 
Method 


1 Introduction 

Free-surface flows of heterogeneous suspensions are abundant in nature and 
technical applications. In principle they are multiphase materials composed of 
a mixture of a liquid and of solid grains of various size. A multitude of in¬ 
teraction mechanisms between these two phases renders the problem of their 
description rather difficult. For example very small grains are bounded to the 
liquid by electrostatic forces, while bigger ones interact mainly by viscous 
forces m- Additionally, inter-grain interactions give rise to the typical com¬ 
plex behavior of granular matter. Often grains have a broad size distribution 
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spanning over several orders of magnitude. Two well known examples are mix¬ 
tures of mud with sand and rocks as well as suspensions of Portland cement, 
sand, and larger aggregates, also known as fresh concrete. While the latter is 
used for construction purposes m, the former gives rise to devastating debris 
flows |23) . 

The simulation of such materials is based either on continuum |351l39j or 
on particle methods [7], depending on whether the investigated effects arise 
from the physics of the fluid or the granular phase. Continuum models are 
appropriate when the rheological behavior of the material can be captured by 
rheometry techniques and phenomenological constitutive laws. However, many 
physical phenomena are eluded by this approach such as size and phase segre¬ 
gation. Examples can be found in concrete casting, where improper mixing or 
vibration leads to inhomogeneities in the physical properties of the hardened 
concrete. In debris flows, size segregation leads to locally changing flow proper¬ 
ties. A flow front rich of large grains with high destructive power is commonly 
observed, followed by a fluid in a more homogeneous tail. Describing this situ¬ 
ation by continuum methods is quite difficult or even physically inappropriate 
[21]. With particle methods, such as the Discrete Element Method (DEM), 
these phenomena can be naturally captured, which makes them an ideal tool 
for the study of the complex behavior of granular materials. 

A complete simulation tool requires a combination of both continuum and 
particle description, which poses serious challenges from a computational point 
of view. If granular and fluid phases are fully coupled, grains represent an ir¬ 
regular and discontinuous boundary for the fluid domain. The relative motion 
of the phases complicates the picture further, because it requires the manage¬ 
ment of continuously evolving interfaces. Eor these reasons, traditional CFD 
solvers such as the standard Finite Element or Finite Volume methods, have 
enormous difficulties to tackle the issue. An attractive alternative is given by 
the Lattice-Boltzmann Method (LBM) |2Ull42j . because of its extreme flexi¬ 
bility in the treatment of elaborate boundary conditions, its ease of imple¬ 
mentation in parallel computing and its superior scaling when compared to 
traditional solvers. For these reasons, much effort has been payed to develop 
of a framework for particle-fluid systems combining the advantages of LBM 
and DEM. The early works in this field are due to Ladd |^[?7l[^ . who first 
coupled LBM and boundaries with imposed velocity. The basic model was en¬ 
hanced by the use of the Immersed Boundary Method [T2]|36l|30l|3T] , by the 
inclusion of turbulence modeling m and extended further to the simulation 
of non-Newtonian rheology models [miiaiin] and for free-surface flows |25l 
HSj . Drawbacks of an approach based on the LBM are its limitation to low- 
Mach and relatively low-Reynolds flows, and the necessity to rely on a regular 
grid, since irregular grids are known to produce a complicated formalism and 
sometimes to lower the accuracy |44j . 

The paper is organized as follows: First a classification of particle suspen¬ 
sions by scales and types of physical interactions is given, before we explain 
how the dynamics of the different phases is addressed. In Sec. [^ we summa¬ 
rize the DEM approach for the granular phase, followed by a section with 
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a comprehensive description of the LBM solver for the fluid phase. Sec. 
explains necessary extensions to the LBM for the simulation of suspensions 
like fluid-particle interaction, non-Newtonian rheology or the representation 
of free surfaces. The experimental validation of the described model completes 
the manuscript in Sec. followed by a brief summary. 


2 Dynamics of Suspensions 

The contribution of grains to the mechanics of the mixture can be of different 
nature depending, among other factors, on the grain size distribution. For a 
phenomenological classification, we use the term small scale when electrostatic 
forces are dominant, medium scale when viscous forces prevail, and large scale 
when inter-particle collisional forces dominate m- Note that the length-scales 
defined by grain size are by no means absolute, but depend on other parameters 
such as the concentration of particles, the viscosity of the liquid, and the state 
of the system, since the same material can exhibit different behaviors when 
sheared at different rates. 

Small scale grain dynamics is governed mainly by interactions of electro¬ 
static nature, e.g. Van der Waals forces. This finer part of the grains, together 
with water, forms a colloidal dispersion. A complete description of this kind of 
material can be found elsewhere |38] . For practical purposes, the mechanics of 
colloidal dispersions is reproduced by continuum methods. A non-Newtonian 
model, however, is generally required, since colloidal dispersions can exhibit 
both shear-thinning behavior and plastic properties. In this work we choose 
to employ the Bingham plastic, a fluid model with a yield stress, well-known 
for its wide applicability [SZIIIT]. It is described as 

J 7 = 0 if fluid does not yield (ct < ay), , , 

\a = ay + Hpij if fluid flows {a > ay), ^ 

where 7 is the magnitude of the shear rate tensor, and Hpi, ay denote plastic 
viscosity and yield stress. In analogy to Newtonian fluids, an apparent viscosity 
(from now on, simply called viscosity) can be locally defined as the ratio of 
shear rate and shear stress 


S'app — 


fJ-pl 


( 2 ) 


As the shear rate 7 approaches zero, the viscosity becomes infinite, giving a 
simple but efficient way to model plastic behavior. 

Medium scale grains are sufficiently big to elude the effects of microscopic 
electrostatic forces and therefore need a different numerical treatment. For 
them the hydrodynamic effects due to the viscous nature of the fluid become 
dominant. In analogy to the smaller scale, grains can be homogenized in the 
fluid. Obtaining an appropriate rheological behavior of the final mixture is 
however more difficult. When experimental data is not available, the value of 
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Grain size [mm] 


Fig. 1 A typical grain size distribution and its effect on the dynamics of the mixture. Big 
grains fall in the collisional regime, small ones in the viscous regime. A transition zone 
exhibits hybrid characteristics. The Bagnold number is calculated with fif = 1.0 Pa s, 
p = 1000 kg/m^, A = 1, 7 = 100 s”! 


the viscosity can be approximated by constitutive relations. A review of these 
models can be found in Ref. El- 

Large scale grain dynamics is dominated by collisions. When collisional 
effects are not damped by viscosity, grains give rise to collective phenomena, 
such as segregation, force percolation or shock waves [Hj. Bagnold defined a 
dimensionless number as the ratio of grain collisional and viscous stresses [ 2 ]. 
It reads 


Ba = 


Psd’^xV^j 

d-f 


(3) 


where Pf is the dynamic viscosity of the liquid, 7 the magnitude of the 
shear rate, ps and ds denote density and characteristic diameter of the grains, 
and Xg their linear concentration (function of the solid fraction Cg as Xg = 
\/[{Cs^max/CgY^^ — 1\ with Cg^max the maximum solid fraction). As illustrated 
in Fig. the Bagnold number is used to distinguish two different regimes, 
where different rheological laws are observed [43]: Mixtures with Ba < 40 are 
dominated by viscosity and therefore the shear stress grows proportionally 
to the shear rate. Mixtures with Ba > 450 are dominated by collisional ef¬ 
fects and grains cannot be homogenized into a continuum description without 
a loss in the descriptive capabilities of the method. An intermediate range 
exists, where both effects are not negligible 

The proposed model follows this classification to efficiently simulate and 
investigate suspensions. Small and medium scale grains are homogenized for a 
fluid formulation with a continuum Bingham model. Large scale grains are rep¬ 
resented by a discrete description. The advantage of this method is that only 
a small portion grains is explicitly represented. This fraction is representative 
both in terms of mass and influence on the rheology of system. 
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3 Dynamics of the granular phase with the DEM 

The granular phase is represented by the DEM, a well-established method for 
granular systems Every grain p is characterized as a Lagrangian element, 
with translation Xp and rotation (ftp as degrees of freedoms. It is subjected 
to multiple interactions that lead to a resultant force Fp and moment Mp. 
These interactions can either be due to collisions, hydrodynamics or volumetric 
forces and are functions of position, orientation and velocity of the particles: 

Fp — Fp ^Xp, Xp, (^p, <^p^ , IVIp — hdp ^Xp, Xp, (^p, <^p^ ■ 

In the simplest case, DEM particles have spherical shape allowing for fast 
contact detection and calculation of the overlap between particles p and 
q namely 



(4) 


Here dp ^ denotes the distance between the center of the spheres with radii 
i?p, Rq. Unfortunately, for most practical applications spheres can only be 
used as a first approximation of the real particle shape. Note that spheres 
exhibit artificial mixing and rolling behavior, which is absent in natural system 
that are not composed of spheres. To overcome these effects we use composite 
elements, created by aggregating a set of spherical particles. While preserving 
the simplicity of the contact calculation, composite elements allow for a more 
realistic representation of granular effects, in particular in the limit of dense 
concentrations. 

Particle-particle interactions are written as the outcome of collisional events 
between particles. Although particles are geometrically described as rigid spheres, 
the overlap ^p^q between particles p and q is used to calculate collisional forces 
and to represent the elastic deformation. We use the law for elastic spheres. 



where Y and v are the Young modulus and the Poisson’s ratio of the mate¬ 
rial, A is a damping constant [5], Reff the effective radius defined as Reff = 
RpRq/{Rp + Rq) and Up ^ the normal vector of the contact surface. The tan¬ 
gential contact force is considered to be proportional to the component of the 
relative velocity of the two spheres laying on the contact surface as 


K,q = i'^rel) - min {vlKelW, fJ-d\\FpJ\) tp^q, 


( 6 ) 


with the tangential shear viscosity coefficient p and the dynamic friction coef¬ 
ficient fJ,d^ thus including Coulomb friction. The tangential unit vector tp ^ is 
obtained normalizing the tangential relative velocity. Wall contacts are calcu¬ 
lated in a similar fashion. 

The time evolution of the system is solved by integrating Newton’s second 
law. 



( 7 ) 
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Fig. 2 a) The regular space discretization of the lattice, b) The 19 discrete velocities allowed 
in the D3Q19 lattice configuration 


While the translational motion is naturally solved in system coordinates, the 
rotational motion requires additional considerations. We use quaternion al¬ 
gebra rather that Euler angles to represent the orientation of the elements, 
and calculate and invert rotation matrices without singularities [8]. Newton’s 
equations in the body-fixed reference frame produce 6 x N scalar equation, 
where N is the number of elements. The Gear predictor-corrector differential 
scheme is used to integrate them M- During the predictor step, tentative val¬ 
ues for particle position, orientation and their derivatives are computed, using 
a Taylor expansion of the previous time step values. The predicted values are 
then used to check for contacts, compute collisional forces and solve Newton’s 
equations of motion. For the corrector step, the difference between the pre¬ 
dicted values for acceleration and their counterpart resulting from Newton’s 
equation is computed. This difference is used to calculate the new corrected 
values for position, orientation and their derivatives. 

The DEM time step should be small enough to resolve the particle contacts. 
If tc is the collision time, then ^ usually < O.ltc- The 

collision time can be estimated for a Hertzian contact as 

= ( 8 ) 

where nieff = mptriq/ {nip + niq) is the effective mass and denotes the 
normal relative velocity at the contact point at the beginning of the collision. 



4 Fluid dynamics with the LBM 

The DEM described in the previous section is coupled with the LBM for solving 
the fluid phase. LBM has been evolving very fast in the last two decades and 
is considered to be one of the most attractive alternatives to traditional CFD 
solvers, especially when problems feature complex boundary conditions. It 
originates from the Boltzmann kinetic theory for the evolution of molecular 
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systems m\- The fluid is described using a distribution function, /(x,c,t), 
defined as the probability density of finding molecules with velocity c at a 
location x and at a given time t. 

In the LBM, the velocity space is discretized by a finite number of ve¬ 
locity vectors, c^, such that /i(x,t) = /(x,Ci,t). We choose to employ the 
D3Q19 lattice cell configuration (3 dimensions and 19 velocities, see Fig. [^, 
which provides the required symmetries to correctly recover the incompress¬ 
ible Navier-Stokes equations. In this work, for simplicity, we use dimensionless 
lattice units {Sx,y,z = 1, d* = 1). 

The reconstruction of macroscopic physical variables such as density pf 
and velocity Uf can be done at every location x and time t by computing the 
first two moments of the distribution function fi (x, t) 

^ * 

The distribution function evolves according to the Lattice Boltzmann Equation 
(LBE), which is written 

fi{x + Ci,t + l) = /*(x,<) -h I7i(x,t). (10) 


Qi represents the collision operator, which in our case corresponds to the linear 
approximation given by Bhatnagar-Gross-Krook [4], 






( 11 ) 


where r is the relaxation time and is the equilibrium distribution function. 
The relaxation time is directly related to the viscosity of the fluid 




r(x,t) - 1/2 
3 


( 12 ) 


For a Newtonian fluid, the relaxation time t is a constant and global param¬ 
eter. However, as stated in Sec. in order to represent most suspensions, a 
non-Newtonian formulation should be employed. To do this, the relaxation 
time is treated as a local variable t(x, t). The equilibrium distribution func¬ 
tion is an expansion in Hermite polynomials of the Maxwell-Boltzmann 
distribution in the limit of small velocities [40]. Using the local macroscopic 
velocity uj and density pf, this yields 


/r(u/.P/) = 


^1 -H3c, - Uf + ^{c,- Uff 



(13) 


where the weights Wi are constants that ensure the recovering of the first and 
second moments of the distribution function (Eq. [^. For the D3Q19 lattice 
configuration they are 

1/3 for f = 1 

1/18 for i = 2...7 

1/36 for f = 8...19. 


Wi = 


(14) 
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To introduce an external force, we employ the scheme developed by Guo 

as 


et al. [19], which consists in modifying Eq. 10 


/* (x + Ci, t + 1) = /, (x, i) + (x, t) + r; (x, t), 


(15) 


where Fi{x,t) is an additional distribution function due to the force field F, 
which can be calculated in a similar fashion as the equilibrium distribution, 

C(x,t) = [3(c, - u/) + 9ci (Cj • Uf)] F. (16) 

With this technique, the computation of the macroscopic velocity field in Eq. 
[^also needs to be modified, 

Uf = ^y]/,(x,t)ci+ F/2^ (17) 

The described approach reproduces the Navier-Stokes equations in the in¬ 
compressible limit. The pressure is directly computable from the density as 

PfifK,t) = pf{x,t)-cl, ( 18 ) 

where Cg is the speed of sound of the fluid, which corresponds to Cg = l/-\/3 
(lattice units). The stability and accuracy of the LBM are guaranteed for small 
Mach numbers, Ma = ||u||/cg <C 1. 

For every time step one first calculates the macroscopic variables, using 
Eq. 1^ and the corresponding equilibrium distribution, from Eq. |13[ Then, one 
uses Eq.[^to evolve the distribution function, which provides the new density 
and velocity of the fluid for the next time step. Being solved mostly at a local 
level, the scheme can be easily implemented in a parallel environment [33) . 


5 Extensions of the LBM for the simulation of suspensions 


To widen the range of applicability of the model to heterogeneous suspensions, 
we need to incorporate a few more features. First, we introduce no-slip moving 


boundaries, necessary for the coupling with the DEM, and second, in Sec. 5.2 


we extend the model to simulate free surfaces. Finally, Sec. |5.3| describes the 
method for non-Newtonian formulations. 


5.1 Coupling with particles 

The coupling with the DEM and the treatment of no-slip boundary conditions 
are performed at a local level by modifying the LBE. Lattice nodes are divided 
into fluid and solid nodes, the latter ones representing particles and walls (see 
Fig.§. Solid nodes are inactive, i.e. on them the LBE is not solved. No-slip 
is performed with the so-called bounce-back rule: every time a distribution 
function fi{x,t) is streaming in the direction i towards a solid node, it gets 
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Fig. 3 Sketch showing how particles or solid objects are discretized on the regular lattice. 
The free-surface is treated in a similar way, with a special type of nodes defining the interface 


reflected back in the opposite direction ik If the boundary is moving, the 
reflected distribution needs to be corrected as 


+ 1) = /*(x,i) - 6w^pu^ ■ Cj, (19) 

where is the local velocity of the wall at the bounce-back location. If the 
wall represents the surface of a particle, the local velocity can be obtained as 

Uu, = Up -I- X ujp, (20) 

where Up and Wp are the linear and angular velocity of the particle, and r^, is 
the vector connecting its center of mass with the bounce-back location. The 
momentum exchange experienced by the reflected distribution can also be used 
to compute the force exerted on the wall when integrated over all bounce-back 
locations. 


Fp = E 


(2/j(x,t) - 6wiPfU^ ■Ci)ci. 


( 21 ) 


Solid boundaries treated this way are located halfway between solid and active 
nodes. This technique was developed for moving boundaries by Ladd [26] and 
Aidun and Lu [T]. 

Particles move over a fixed, regular grid. Of course the node classification 
into fluid and solid is not fixed but needs to be updated. Following the particle 
motion, fluid nodes are created (deleted) in the wake (front) of moving parti¬ 
cles. The macroscopic density and the velocity of the newly created nodes are 
calculated as the average over the values in the neighborhood as initial values 
for the distribution function /i(x, t) through Eq. 13 Deleted fluid nodes are 


converted to solid ones and therefore made inactive. Both processes introduce 
small variations in the global mass and momentum. However, due to the fact 
that all our simulations are performed in the incompressible limit (variation of 
density are very small), and that fluid nodes close to a particle possess nearly 
the same velocity as the particle, we expect these variations to be negligible. 
Another problem is the representation of the particle boundaries on the reg¬ 
ular lattice, which leads to a zig-zag approximation of the spherical shapes. 
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An alternative way to overcome these problems is the use of the Immersed 
Boundary Method [12] or of a fictitious domain nmn]. Both methods are 
more precise and smooth the ill effects of particles traveling though the lat¬ 
tice. At the same time, they require additional computations and are therefore 
avoided following the spirit of this paper. 

When two particles approach each other, the distance between the surfaces 
can become smaller than the lattice node spacing, resulting in an imprecise 
resolution of the collision process. To overcame this problem, we use the lu¬ 
brication theory of Nguyen and Ladd |34j . In this theory, when two particle 
are moving with a relative velocity Urei, the correction force 



( 22 ) 


is added, where Sp^q = —^p,q is the distance between the particle surfaces and 
diub denotes a cut-off distance above which no force is computed. 

5.2 Free surface representation 

We employ the mass tracking algorithm described in Refs. wm which, despite 
its simplicity, leads to a stable and accurate surface evolution. Fluid nodes are 
further divided into liquid, interface and gas nodes: Liquid and interface nodes 
are considered active, and the LBE is solved. The remaining nodes are the gas 
nodes and are inactive, with no evolution equation. Liquid and gas nodes are 
never directly connected, but through an interface node (see Fig. |^. 

An additional macroscopic variable for the mass m/ (x, t) stored in a node 
is required, dehned as 



(23) 


The mass is updated using the equation 


mf{x,t + l) = mf{x,t) [fi>{x + Ci,t) - /*(x,t)], (24) 


where ai is a parameter determined by the nature of the neighbor node in the 
i direction, 

{ i [m/(x, t) +mf{'x + Ci,t)] if the neighbor node is interface, 

1 if the neighbor node is liquid, (25) 

0 if the neighbor node is gas. 

When the mass becomes zero (m/(x, t) = 0), the interface node is transformed 
into gas, with all liquid nodes connected to it becoming interface. Analogously, 
an interface node whose mass reaches the density (m/(x, t) = Pf{x, t)) is trans¬ 
formed into liquid, and all connected gas nodes become interface. However, 
due to the discrete integration, these equalities are not in general satisfied. 
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Fig. 4 Representation of the rheology model employed for plastic fluids. The approximation 
of the Bingham model is limited by the maximum and minimum acceptable values for the 
relaxation time r and therefore for the viscosity /ry 


The surplus of mass is equally distributed to the neighboring interface nodes, 
conserving the total mass of the system. 

Because gas nodes are not active, there are no distribution functions stream¬ 
ing from gas nodes to interface nodes. These missing distribution functions are 
computed from the macroscopic variables at the interface, atmospheric density 
Patm and interface velocity Uint, as 

/p (x-b Cp , t -h 1) = /f'^(u„t,Patm) + fP{Viint,Patm) “ (26) 

Note that this implies that gas nodes have the same macroscopic velocity as 
the connected interface nodes. 


5.3 Bingham plastic rheology model 

The presence of the small particle fraction in the fluid leads to non-Newtonian 
behavior, that needs to be considered. For the LBM this implies that the 
relaxation time r is not a global parameter for the system, but rather r = 
t(x, t). A non-linear dependency of viscosity, and thus of r, on the shear rate 
requires an explicit computation of the shear rate tensor. This can be done with 
ease in the LBM from the non-equilibrium part of the distribution functions, 

7a&(x, t) = ^haCi,b (/*(X, t) - /f‘*(x, t)) . (27) 

With the second invariant of the shear rate tensor 

A ^) = X! X! AAab, (28) 

a b 

the magnitude of the shear rate is calculated as 


7(x,t) = y^2r^(x,t) 


(29) 
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Fig. 5 Rheology test on the fresh cement paste. A linear Bingham approximation is used 
to fit the data, obtaining fipi and cry 


This can be included in any constitutive equation for purely viscous fluids. As 
outlined before (Sec.[^, we choose the Bingham constitutive model and get a 
new form of Eq. [^for the explicit update of t, 

r(x,f) = y3(p,, + ^), (30) 

The accuracy and stability of LBM are guaranteed only over a certain range 
of values for t. This limits the applicability of Eq. because r diverges 
when 7 —0. Following Svec et al. [15], we use a simple solution to this 
problem, imposing that < r^ax- Reasonable values for Tmin 

and Tmax are, respectively, 0.501 and 3.5. The constitutive equation arising 
from this approach is that of a tri-viscosity fluid (see Fig. |^. If < fJ-pi 

the model represents a bi-viscosity, and if /r f^max , the approximation of 

the Bingham model is fair. With these extensions the model is complete and 
we can address examples. 


6 Experimental validation by a gravity-driven flow 

The capabilities of the model are shown by comparing with an experiment, 
featuring a free-surface flow of a suspension under the effect of gravity. We 
employ fresh concrete, since it poses all the challenges necessary to validate 
the method: a non-Newtonian rheology and an irregular granular phase. The 
cement paste is obtained with a commercial Portland cement of type CEM 
I 42.5N. Water is added until a water/cement ratio of 0.4 is reached. The 
rheology of the obtained paste is measured with a coaxial rotational viscometer 
Haake RV20. The measurement procedure consists in the uniform shearing of 
the paste at 200 s“^ for 120 s, followed by a shear rate continuous ramp from 
0 up to 200 s“^ occurring over 120 s [Hj. The obtained rheological curve is 
shown in Fig. 
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Fig. 6 Experiment setup with an inclined plane of 933 X 700mm. The internal size of the 
box is 150 X 150 mm. The dark gray area represents the surface covered by the sandpaper 


The paste is mixed with 1000 silica rounded pebbles with radius R = 
4.0 -h 8.0 mm. The total weight of the grains is 2.629 kg, and the density 
Ps = 2680 kg/m^. The components are mixed in a bowl until homogenization 
and then vibrated for degassing. The final mixture is poured in a 150 x 150 mm 
rectangular box, open on top and bottom and positioned over a wooden board 
inclined at 15°. The board surface is upholstered with sandpaper and wetted 
before the start of the test. The test is performed by steadily lifting the box, 
and letting the sample spread on the board under the sole effect of gravity. The 
flow falls in the intermediate regime (see Fig.[^. Collisional effects are therefore 
not dominant, but still important. The geometry of the test is illustrated in 
Fig.§ and Fig. (a) is a picture of the final deposition of the sample. 

The same environment is set up with a simulation on an LBM lattice of 
350 X 250 X 80 nodes, with the lattice spacing corresponding to 2.0 • 10“^m 
in physical units. The initial configuration of the fluid is a cube with edge 
length of 0.15m, corresponding to 75 x 75 x 75 liquid nodes. The pebbles 
are represented with 1000 discrete elements, each composed of 4 spheres with 
tetrahedral structure. The total number of spheres is 4000. The box is repre¬ 
sented by a set of moving walls and is set as solid boundary both for fluid and 
granular solvers. The lifting speed of the box is 0.15m/s. The properties of the 
fluid are obtained from the viscometer data, as represented in Fig. A good 
fit is obtained with a Bingham model with plastic viscosity Ppi = 0.15 Pa • s 
and yield stress ay = 62 Pa. The model is imprecise for lower shear rates, 
which is one of the limitations of the chosen linear approach. Fig. shows 
the results of the simulation on the longitudinal cross section of the sample. 
The evolution of the shear rate and the particle distribution can be tracked 
continuously. The final shapes of the experimental and numerical solution are 
compared in Fig.[7](c), showing excellent agreement. 

A good compromise between stability and speed is obtained with a time 
step of = 3.0 • 10“® s. This sets the maximum allowable speed in 
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Fig. 7 Final shape of the flowing mass. a,d) Numerical shape; b,e) Experimental shape. 
Fluid mass opacity in the numerical shape is lowered for a better visualization of particles; 
c,f) Comparison of numerical shape (solid line) and experimental shape (dashed line). The 
background grid has 5cm spacing 


the system as 0.667 m/s. The parameters in lattice units are then viscos- 
ity = 6.25 • 10“"^, relaxation time = 7.75 • 10“®, and gravity 

||gLBM|| _ ^ ■ 10“®. The simulation is stopped when 95% of the fluid has 

reached the maximum viscosity. The total simulation time is 63 hours, with a 
parallel run on 4 cores with an Intel Xeon E5-1620 processor at 3.60 GHz.. 


7 Summary 

In this paper a model for the simulation of the flow of suspensions was pro¬ 
posed. The multiscale nature of the model is justified by the different interac¬ 
tion mechanisms acting between the liquid and the granular phase. A practical 
mean of phenomenological classification of interactions is given by the Bagnold 
number: Small grains are considered to be governed by the viscous nature of 
the liquid and are modeled as part of the fluid phase itself with the use of 
a plastic non-Newtonian formulation. Grains with a sufficiently large size are 
dominated by collisional mechanisms. This is modeled with a two-way coupling 
between fluid and grains, along with the resolution of particle contacts. 
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Fig. 8 Dynamic viscosity contour on the longitudinal cross section of the simulation. Par¬ 
ticles are represented in light gray and walls in dark gray. The yielded region of the fluid 
grows from a small portion close to the box walls to the whole sample during the first part 
of the simulation. When a new equilibrium is reached, the yielded region reduces and the 
flow is slowed 


The problem was solved with a hybrid of the Discrete Element Method 
for grains and the Lattice-Boltzmann Method for fluids. A combination of the 
most successful advances in these methods was employed. The mass-tracking 
algorithm allows an inexpensive way to simulate free surfaces, while the vari¬ 
able relaxation time formulation can reproduce non-Newtonian constitutive 
laws. The hydrodynamic interactions with the granular phase were fully solved 
with the bounce-back rule for coupling non-slip moving boundaries and fluid. 
The proposed model finds its best application in the simulation of real flows 
and in particular of heterogeneous suspensions with a granular phase that fea¬ 
tures a complete size distribution, due to its multiscale nature. The intrinsic 
advantages of the Lattice-Boltzmann solver, with its high- level performance 
and its relatively simple implementation make it a good choice for the fast 
development of such methods. Moreover, the core of the solver works at a 
local level, making the parallelization of the code easy and natural. Grain- 
grain interactions were solved with a Discrete Element Method. We assured 
that the scaling of the particle solver was not too far from the almost linear 
performances of the fluid solver. The Hertzian contact law was used, and a 
formulation for non-spherical particles was included. The capabilities of the 
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approach were shown by comparing to an experimental free-surface flow of a 
fresh concrete sample. An excellent agreement between numerical and experi¬ 
mental data was found in the comparison of the final shape of the sample. The 
results of the simulation can provide insight into the mechanics of the flow. 
The spatial distribution of particles can be tracked, along with the variables 
of the flow: velocity, pressure, shear rate and viscosity. 

Another challenging application of the model is the prediction of debris 
flows, which can hardly be assessed experimentally. Future works will focus 
on the rheology of debris materials and on the full simulation of events for 
deeper physical understanding, and on techniques for the design of effective 
protection measures. 
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